﻿#pragma once

#include <functional>

#include <xy/alg/interpolator.h>
#include <xy/math/tensor.h>

namespace xy
{

namespace alg
{

class integrator
{
  public:
    /**
     * @brief 构造函数
     *
     * @param[in] N 节点数
     */
    explicit integrator(std::size_t N);

    /**
     * @brief 获取节点数
     *
     * @return std::size_t
     */
    std::size_t node_num() const;

    /**
     * @brief 积分接口
     *
     * @param[in] f
     * @param[in] a
     * @param[in] b
     * @return double
     */
    virtual double solve(std::function<double(double)> f, double a, double b) = 0;

  protected:
    std::size_t m_N; // 节点数
};

struct integrator_trapezoidal : public integrator
{
    /**
     * @brief 节点数
     *
     * @param[in] N
     */
    explicit integrator_trapezoidal(std::size_t N);

    /**
     * @brief 积分接口
     *
     * @param[in] f
     * @param[in] a
     * @param[in] b
     * @return double
     */
    double solve(std::function<double(double)> f, double a, double b) override;
};

struct integrator_simpson : public integrator
{
    /**
     * @brief 节点数，必须为奇数
     *
     * @param[in] N
     */
    explicit integrator_simpson(std::size_t N);

    /**
     * @brief 积分接口
     *
     * @param[in] f
     * @param[in] a
     * @param[in] b
     * @return double
     */
    double solve(std::function<double(double)> f, double a, double b) override;
};

struct integrator_gauss : public integrator
{
    /**
     * @brief 节点数
     *
     * @param[in] N
     */
    explicit integrator_gauss(std::size_t N);

    /**
     * @brief 积分接口
     *
     * @param[in] f
     * @param[in] a
     * @param[in] b
     * @return double
     */
    double solve(std::function<double(double)> f, double a, double b) override;

    std::vector<double> weights; // 积分权
    std::vector<double> nodes;   // 积分节点
};

template <std::size_t Dim> struct integrator_newton_tensor
{
    using Grid = math::tensor<Dim, double>;
    using Param = geom::vec<Dim, double>;
    math::tensor<Dim, double> weights; // 积分权

    /**
     * @brief Construct a new integrator newton tensor object
     *
     * @param[in] N
     */
    explicit integrator_newton_tensor(std::size_t N) : weights(N + 1, 0)
    {
        // 计算权的分布
        math::tensor<1, double> baseWeights(N + 1, 0);

        // 使用插值多项式积分计算权
        std::vector<double> x(N + 1);
        for (std::size_t i = 0; i < N + 1; i++)
            x[i] = i;

        for (std::size_t i = 0; i < N + 1; i++)
        {
            std::vector<double> f(N + 1, 0);
            f[i] = 1;

            interpolator_newton<double> interpolator(x, f);
            baseWeights.data()[i] = interpolator.interpolate().integrate(0, N);
        }
        weights = power<Dim>(baseWeights);
    }

    /**
     * @brief 积分接口，计算函数 f: R^Dim -> Any 在区域 [lower, upper] 上的积分值
     *
     * @param[in] f
     * @param[in] lower
     * @param[in] upper
     * @return vecxd
     */
    math::vecxd solve(std::function<math::vecxd(const Param &)> f, const Param &lower, const Param &upper)
    {
        // 离散 Dim 维区域为网格，并计算单元体积
        math::tensor<Dim, double> grid(weights.base(), 0);
        double volume = 1;
        for (std::size_t i = 0; i < Dim; i++)
            volume *= (upper[i] - lower[i]) / (grid.base() - 1);

        // 获得输出向量的长度
        std::size_t n = f(lower).size();
        math::vecxd I(n);

        // 由于输出分量之间无关，因此可以并行计算
        for (std::size_t i = 0; i < n; i++)
        {
            do
            {
                Param param;
                auto base = grid.base();
                for (std::size_t i = 0; i < Dim; ++i)
                {
                    double h = (upper[i] - lower[i]) / (base - 1);
                    param[i] = lower[i] + grid.index()[i] * h;
                }

                // 对每个网格点计算函数值并保存
                auto y = f(param);
                grid.at_index() = y[i];
            } while (grid.move_index());
            I[i] = weights * grid * volume;
        }
        return I;
    }
};

} // namespace alg

} // namespace xy
